1 Data import

The data is included within this reproducible report and can be downloaded here (object 'rbt_public_1') and here (object 'rbt_public_0').

1.1 Data wrangling

# wrangling
rbt_public_l <- pivot_longer(rbt_public_1, 
                             cols = (abs1_tsm_1:abs2_tch_5), 
                             names_to = "variable", 
                             values_to = "value") %>%
    mutate(treat = case_when(                    # create treatment variable
        str_detect(variable, "abs1_") ~ toupper(treat1),
        str_detect(variable, "abs2_") ~ toupper(treat2)),
        variable = str_sub(variable, 6, -1)) # delete title page substring

# pivot them into wide format
rbt_public_w <- pivot_wider(rbt_public_l,
                            names_from = variable,
                            values_from = value,
                            id_cols = c(session, treat))

# invert trust items & build scales
inv7fct <- function(x) (8-as.numeric(x))
data_rbt <- rbt_public_w %>%
    mutate_at(vars(tru_exp_1:tru_ben_4),
              list(~inv7fct(.))) %>% # recoding 1-->7, 2-->6, ...
    rename(exp_1 = tru_exp_1, 
           exp_2 = tru_exp_2, 
           exp_3 = tru_exp_3, 
           exp_4 = tru_exp_4, 
           exp_5 = tru_exp_5, 
           exp_6 = tru_exp_6,
           int_1 = tru_int_1, 
           int_2 = tru_int_2, 
           int_3 = tru_int_3, 
           int_4 = tru_int_4,
           ben_1 = tru_ben_1, 
           ben_2 = tru_ben_2, 
           ben_3 = tru_ben_3, 
           ben_4 = tru_ben_4) %>% 
    group_by(session) %>% 
    mutate(meas_rep = 1:n()) %>% 
    ungroup() %>% 
    full_join(., rbt_public_0)
## Joining, by = "session"
# data frames with sum scales
data_scales <- data_rbt%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
         TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))

data_scales_lables <- data_rbt%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Treatment = fct_recode(Treatment,
                                `Colored badges` = "CB",
                                `Control Condition` = "CC",
                                `Greyed out badges` = "GB"),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
        `Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), 
                                               na.rm = T))

data_scales_wide <- data_scales%>%
  select(Experitse, Integrity, Benevolence,
         TSM, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

data_scales_lables_wide <- data_scales_lables%>%
  select(Experitse, Integrity, Benevolence,
         `Topic specific multiplism`, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

# Data of treatment check
data_tch <- data_rbt %>% 
  select(starts_with("tch"), meas_rep, treat)

data_tch_n <- data_tch%>%
  mutate_all(function(x) ifelse(x == -999, NA, x)) %>% 
  filter(rowSums(is.na(data.frame(tch_1, tch_2, tch_3, tch_4))) != 4)

PID_list <- data_rbt$session %>% unique()

2 Sample

2.1 Sample size

length(PID_list)
## [1] 257
data_rbt %>%
  filter(meas_rep == 2) %>%
  select(session) %>% 
  distinct() %>% 
  nrow()
## [1] 257

2.2 Skipped after first measurement

sum(is.na(data_rbt%>%
            filter(meas_rep == 2)%>%
            pull(exp_1)))
## [1] 0

2.3 Demographics

data_rbt %>% 
  filter(meas_rep == 2)%>%
  mutate(education = as.factor(education),
         sex = as.factor(sex),
         age = as.factor(age)) %>% 
  select(age, education, sex) %>% 
  skim(.)
Data summary
Name Piped data
Number of rows 257
Number of columns 3
_______________________
Column type frequency:
factor 3
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
age 0 1 FALSE 3 50: 108, 16: 78, 35: 71
education 0 1 FALSE 13 app: 61, L34: 56, L12: 26, L12: 19
sex 0 1 FALSE 2 f: 134, m: 123

2.3.1 Count on sex

data_rbt %>%
  filter(meas_rep == 2)%>%
  pull(sex) %>% 
  table(.)
## .
##   f   m 
## 134 123

2.3.2 Count on education

data_rbt %>%
  filter(meas_rep == 2)%>%
  pull(education) %>% 
  table(.)
## .
## app_11 app_12 app_13  app_5  L12_1  L12_2  L12_3  L12_4 L34_10  L34_6  L34_7 
##     11      3     61      6     19      9     26     18     14      8     18 
##  L34_8  L34_9 
##     56      8

2.3.3 Count on age

data_rbt %>%
  filter(meas_rep == 2)%>%
  pull(age) %>% 
  table(.)
## .
##  16  35  50 
##  78  71 108

3 Psychometric properties of the measurements

3.1 Descriptive overview

skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))

my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))

data_scales%>%
  my_skim(.)
Data summary
Name Piped data
Number of rows 514
Number of columns 37
_______________________
Column type frequency:
character 4
factor 1
numeric 32
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
session 0 1 64 64 0 257 0
treat 0 1 2 2 0 3 0
sex 0 1 1 1 0 2 0
education 0 1 5 6 0 13 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Treatment 0 1 FALSE 3 GB: 183, CC: 170, CB: 161

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist skewn kurto maxabszvalue
tsm_1 0 1 2.52 0.79 1 2.00 3.00 3.00 4 ▂▇▁▇▂ -0.11 -0.43 1.92
tsm_2 0 1 2.72 0.82 1 2.00 3.00 3.00 4 ▁▇▁▇▃ 0.05 -0.74 2.11
tsm_3 0 1 2.56 0.83 1 2.00 3.00 3.00 4 ▂▇▁▇▂ 0.02 -0.58 1.88
tsm_4 0 1 2.82 0.82 1 2.00 3.00 3.00 4 ▁▅▁▇▃ -0.21 -0.56 2.22
tsc_1 0 1 2.65 0.79 1 2.00 3.00 3.00 4 ▁▆▁▇▂ -0.16 -0.40 2.08
tsc_2 0 1 2.30 0.81 1 2.00 2.00 3.00 4 ▂▇▁▆▁ 0.19 -0.45 2.11
tsc_3 0 1 2.67 0.80 1 2.00 3.00 3.00 4 ▁▆▁▇▂ -0.20 -0.39 2.10
exp_1 0 1 5.08 1.42 1 4.00 5.00 6.00 7 ▁▁▅▅▇ -0.54 -0.16 2.87
exp_2 0 1 5.16 1.41 1 4.00 5.00 6.00 7 ▁▁▅▃▇ -0.45 -0.44 2.95
exp_3 0 1 5.16 1.42 1 4.00 5.00 6.00 7 ▁▁▅▃▇ -0.53 -0.19 2.92
exp_4 0 1 5.15 1.46 1 4.00 5.00 6.00 7 ▁▂▃▃▇ -0.53 -0.35 2.84
exp_5 0 1 5.04 1.38 1 4.00 5.00 6.00 7 ▁▁▅▅▇ -0.40 -0.28 2.93
exp_6 0 1 5.01 1.46 1 4.00 5.00 6.00 7 ▁▂▅▃▇ -0.42 -0.47 2.75
int_1 0 1 4.94 1.42 1 4.00 5.00 6.00 7 ▁▂▆▅▇ -0.38 -0.34 2.76
int_2 0 1 4.99 1.39 1 4.00 5.00 6.00 7 ▁▂▆▅▇ -0.29 -0.50 2.87
int_3 0 1 4.73 1.33 1 4.00 5.00 6.00 7 ▁▂▇▅▇ -0.24 0.05 2.80
int_4 0 1 4.89 1.34 1 4.00 5.00 6.00 7 ▁▂▇▆▇ -0.29 -0.16 2.90
ben_1 0 1 4.80 1.36 1 4.00 5.00 6.00 7 ▁▂▇▅▇ -0.21 -0.25 2.80
ben_2 0 1 4.87 1.38 1 4.00 5.00 6.00 7 ▁▂▇▅▇ -0.28 -0.25 2.80
ben_3 0 1 4.89 1.42 1 4.00 5.00 6.00 7 ▂▂▆▅▇ -0.35 -0.39 2.74
ben_4 0 1 4.85 1.32 1 4.00 5.00 6.00 7 ▁▂▆▇▇ -0.24 -0.28 2.92
tch_1 0 1 -176.59 384.37 -999 1.00 3.00 3.00 4 ▂▁▁▁▇ -1.67 0.79 2.14
tch_2 0 1 -197.98 401.39 -999 1.00 3.00 4.00 4 ▂▁▁▁▇ -1.49 0.23 2.00
tch_3 0 1 -264.33 443.31 -999 -999.00 2.00 3.00 4 ▃▁▁▁▇ -1.05 -0.89 1.66
tch_4 0 1 -159.05 368.96 -999 1.00 3.00 3.00 4 ▂▁▁▁▇ -1.83 1.37 2.28
tch_5 0 1 -240.87 430.18 -999 1.00 2.00 3.00 4 ▂▁▁▁▇ -1.19 -0.58 1.76
meas_rep 0 1 1.50 0.50 1 1.00 1.50 2.00 2 ▇▁▁▁▇ 0.00 -2.00 1.00
age 0 1 35.54 14.29 16 16.00 35.00 50.00 50 ▆▁▅▁▇ -0.34 -1.50 1.37
Experitse 0 1 5.10 1.27 1 4.17 5.00 6.00 7 ▁▁▆▆▇ -0.42 -0.15 3.22
Integrity 0 1 4.88 1.19 1 4.00 4.75 5.75 7 ▁▂▇▇▅ -0.18 -0.02 3.26
Benevolence 0 1 4.85 1.20 1 4.00 4.75 5.75 7 ▁▂▇▇▅ -0.13 -0.28 3.21
TSM 0 1 2.65 0.57 1 2.25 2.75 3.00 4 ▁▂▇▃▂ 0.08 0.16 2.92

3.2 METI

3.2.1 Dimensionality (CFA)

First we specified two consecutive threedimensional CFA models

# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                             int_1 + int_2 + int_3 + int_4 +
                             ben_1 + ben_2 + ben_3 + ben_4"

cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 2))

cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4"

cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4"

cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_rbt%>%filter(meas_rep == 1))
## Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
##                 is not positive definite;
##                 use lavInspect(fit, "cov.lv") to investigate.
# define a function which prints the fit
fpf <- function(x){  
  fm_tmp <- fitmeasures(x)
  return(cat(sprintf(
          "χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
           round(fm_tmp[c("chisq")],3), 
                 fm_tmp[c("df")],
           round(fm_tmp[c("cfi")],3),
           round(fm_tmp[c("tli")],3),
           round(fm_tmp[c("rmsea")],3),
           round(fm_tmp[c("srmr")],3),
           round(fm_tmp[c("srmr_between")],3),
           round(fm_tmp[c("srmr_within")],3))))
}

# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)

χ2 = 308.08, df = 77, CFI = 0.931, TLI = 0.918, RMSEA = 0.108, SRMR = 0.04, SRMRbetween = NA, SRMRwithin = NA

fpf(cfa_meti_1)

χ2 = 194.143, df = 74, CFI = 0.964, TLI = 0.956, RMSEA = 0.079, SRMR = 0.031, SRMRbetween = NA, SRMRwithin = NA

cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)

χ2 = 297.077, df = 77, CFI = 0.941, TLI = 0.93, RMSEA = 0.105, SRMR = 0.033, SRMRbetween = NA, SRMRwithin = NA

fpf(cfa_meti_2)

χ2 = 196.021, df = 74, CFI = 0.967, TLI = 0.96, RMSEA = 0.08, SRMR = 0.025, SRMRbetween = NA, SRMRwithin = NA

anova(cfa1d_meti_1, cfa_meti_1)
## Chi-Squared Difference Test
## 
##              Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_1   74 9528.6 9638.6 194.14                                  
## cfa1d_meti_1 77 9636.5 9735.9 308.08     113.94       3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cfa1d_meti_2, cfa_meti_2)
## Chi-Squared Difference Test
## 
##              Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_2   74 8883.3 8993.3 196.02                                  
## cfa1d_meti_2 77 8978.4 9077.8 297.08     101.06       3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In an next step we ran a two-level CFA …

# onedimensional model
mcfa1d_meti_model <- "level: 1
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4
                    
                    level: 2
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4"



mcfa_meti_model <- "level: 1
                    exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_w =~ int_1 + int_2 + int_3 + int_4 
                    ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
                    
                    level: 2
                    exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_b =~ int_1 + int_2 + int_3 + int_4 
                    ben_b =~ ben_1 + ben_2 + ben_3 + ben_4"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_rbt, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_rbt, cluster = "session")
fpf(mcfa1d_meti)

χ2 = 437.579, df = 154, CFI = 0.953, TLI = 0.944, RMSEA = 0.06, SRMR = 0.159, SRMRbetween = 0.052, SRMRwithin = 0.107

fpf(mcfa_meti)

χ2 = 313.519, df = 148, CFI = 0.972, TLI = 0.966, RMSEA = 0.047, SRMR = 0.089, SRMRbetween = 0.04, SRMRwithin = 0.049

anova(mcfa1d_meti, mcfa_meti)
## Chi-Squared Difference Test
## 
##              Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## mcfa_meti   148 18217 18539 313.52                                  
## mcfa1d_meti 154 18329 18626 437.58     124.06       6  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2.2 Reliability (McDonalds \(\omega\))

# First Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("exp")))$est
## [1] 0.9492506
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("int")))$est
## [1] 0.8795378
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("ben")))$est
## [1] 0.8819123
# Second Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("exp")))$est
## [1] 0.9504771
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("int")))$est
## [1] 0.900965
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("ben")))$est
## [1] 0.9137172

3.3 Topic specific multiplism

3.3.1 Dimensionality (CFA)

First we specified two consecutive onedimensional CFA models

cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
                  tsm_2 ~~ tsm_4"

cfa_tsm_1 <-
  cfa(cfa_tsm_model, 
      data = data_rbt %>% 
        filter(meas_rep == 1),
      missing = "fiml")
fpf(cfa_tsm_1)

χ2 = 3.991, df = 4, CFI = 1, TLI = 1, RMSEA = 0, SRMR = 0.03, SRMRbetween = NA, SRMRwithin = NA

cfa_tsm_2 <-
  cfa(cfa_tsm_model, 
      data = data_rbt %>% 
        filter(meas_rep == 2))
fpf(cfa_tsm_2)

χ2 = 5.472, df = 4, CFI = 0.988, TLI = 0.982, RMSEA = 0.038, SRMR = 0.05, SRMRbetween = NA, SRMRwithin = NA

In an next step, we ran a two-level CFA …

mcfa_tsm_model <- "level: 1
                    tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ tsm_4
                    

                    level: 2
                    tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ 0*tsm_2
tsm_1 ~~ tsm_2"

mcfa_tsm <- cfa(mcfa_tsm_model, data = data_rbt, cluster = "session")
fpf(mcfa_tsm)

χ2 = 4.328, df = 3, CFI = 0.995, TLI = 0.98, RMSEA = 0.029, SRMR = 0.092, SRMRbetween = 0.071, SRMRwithin = 0.02

3.3.2 Reliability (McDonalds \(\omega\))

# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==1) %>% select(starts_with("tsm")))$est
## [1] 0.6962778
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==2) %>% select(starts_with("tsm")))$est
## [1] 0.6046964

3.4 Treatment check

3.4.1 Dimensionality (CFA)

We specified two consecutive onedimensional CFA models

cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5
                  tch_1 ~~ tch_2"

cfa_tch_model2 <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5"

cfa_tch_1 <- cfa(cfa_tch_model, 
                 data = data_tch_n%>%
                   filter(meas_rep == 1))
fpf(cfa_tch_1)

χ2 = 3.982, df = 4, CFI = 1, TLI = 1, RMSEA = 0, SRMR = 0.012, SRMRbetween = NA, SRMRwithin = NA

cfa_tch_2 <- cfa(cfa_tch_model2, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)

χ2 = 8.846, df = 5, CFI = 0.996, TLI = 0.991, RMSEA = 0.066, SRMR = 0.016, SRMRbetween = NA, SRMRwithin = NA

3.4.2 Reliability (Ordinal McDonalds \(\omega\))

# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 1) %>% select(starts_with("tch")))$est
## [1] 0.8354668
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est
## [1] 0.9357098

3.5 Table of (M)CFA fit-indices

tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
  column_to_rownames(var = "rownames")%>%
  knitr::kable(., digits = 3)
1d CFA METI 1 1d CFA METI 2 3d CFA METI 1 3d CFA METI 2 1d MCFA METI 3d MCFA METI 1d CFA TSM 1 1d CFA TSM 2 1d MCFA TSM 1d CFA TCH 1 1d CFA TCH 2
χ2 308.080 297.077 194.143 196.021 437.579 313.519 3.991 5.472 4.328 3.982 8.846
df 77.000 77.000 74.000 74.000 154.000 148.000 4.000 4.000 3.000 4.000 5.000
CFI 0.931 0.941 0.964 0.967 0.953 0.972 1.000 0.988 0.995 1.000 0.996
TLI 0.918 0.930 0.956 0.960 0.944 0.966 1.000 0.982 0.980 1.000 0.991
RMSEA 0.108 0.105 0.079 0.080 0.060 0.047 0.000 0.038 0.029 0.000 0.066
SRMR 0.040 0.033 0.031 0.025 0.159 0.089 0.030 0.050 0.092 0.012 0.016
SRMRbetween NA NA NA NA 0.052 0.040 NA NA 0.071 NA NA
SRMRwithin NA NA NA NA 0.107 0.049 NA NA 0.020 NA NA
BIC 9735.922 9077.757 9638.632 8993.349 18625.957 18539.351 2421.423 2341.021 4734.532 1709.467 1853.905
AIC 9636.548 8978.383 9528.610 8883.327 18329.002 18216.942 2385.932 2319.726 4645.445 1676.277 1822.088

4 Results of the treatmentcheck

4.1 Plot

res_tch_data <- data_tch%>%
  gather(variable, value, starts_with("tch_"))%>%
  group_by(treat, variable, value)%>%
  mutate(value = as.character(value)) %>% 
  summarize(freq = n())%>%
  ungroup()%>%
  mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
                           treat == "CB" ~ "Colored badges",
                           treat == "CC" ~ "Control condition",
                           T ~ treat),
         value = case_when(value =="-999" ~ "don't know",
                           T ~ value),
         variable = case_when(variable == "tch_1" ~ "item 1", 
                              variable == "tch_2" ~ "item 2", 
                              variable == "tch_3" ~ "item 3", 
                              variable == "tch_4" ~ "item 4", 
                              variable == "tch_5" ~ "item 5"),
         Frequency = freq)
## `summarise()` regrouping output by 'treat', 'variable' (override with `.groups` argument)
res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) + 
  geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
  scale_size_continuous(range = c(3,15)) + 
  scale_color_gradient(low = "grey95", high = "grey65") +
  guides(color=guide_legend(), size = guide_legend()) +
  facet_wrap(~treat) + 
  theme_ipsum_ps() + 
  ggtitle("Results of the treatment check", "Frequency per item and experimental condition") + 
  ylab("") + 
  xlab("")
  
res_tch_plot

res_tch_plot_pub <- ggplot(res_tch_data %>%
                         mutate(treat = case_when(treat == "Greyed out badges" ~ "Grey badges",
                                                  TRUE ~ treat)), aes(variable, value)) + 
  geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
  scale_size_continuous(range = c(3,15)) + 
  scale_color_gradient(low = "grey95", high = "grey65") +
  guides(color=guide_legend(), size = guide_legend()) +
  facet_wrap(~treat) + 
  theme_ipsum_ps() + 
  ylab("") + 
  xlab("")

#ggsave("Fig/res_tch_plot_public.svg", res_tch_plot, width = 11, height = 6)
#ggsave("8_reproducible_documentation_of_analyses/Fig/Fig7.jpg", res_tch_plot_pub, width = 130, height = 50, units = "mm", dpi = 300, scale = 2.4)

4.2 Effect size

res_tch_data_A <- data_tch_n%>%
  filter(treat != "CC")%>%
  filter(tch_1 != -999)

effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)   
## 
## Vargha and Delaney A
## 
## A estimate: 0.72609 (medium)

5 Graphical exploration

5.1 Plot Hyp 1

data_scales_lables%>%
  gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
  ggplot(., aes(x = Treatment, y = Value)) + 
  geom_violin(adjust = 1.5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  facet_wrap(~ Variable, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 1)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  ylim(1,7) +
  hrbrthemes::theme_ipsum_ps()

# Export data for combined plot
#data_scales_lables %>% 
#  mutate(Study = "Study 3") %>% 
#  write_csv(., "8_reproducible_documentation_of_analyses/Fig/data_scales_lables_study_3.csv")

5.2 Descriptive Effect Sizes Hyp 1

A_GB_CC <- data_scales_lables%>%
  filter(treat != "CB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_CC_CB <- data_scales_lables%>%
  filter(treat != "GB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_GB_CB <- data_scales_lables%>%
  filter(treat != "CC")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate
GB CC
CC A = 0.56, d = 0.21
CB A = 0.56, d = 0.2 A = 0.49, d = -0.02

5.3 Hyp 2/3

plot_hyp2_1 <- ggplot(data_scales_lables, 
                      aes(x=`Topic specific multiplism`, y = Integrity)) + 
  geom_jitter() +
  facet_wrap(~ Treatment, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 2/3)",
       subtitle = "Jitter plot per treatment") +
  hrbrthemes::theme_ipsum_ps()


plot_hyp2_1 + stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plot_hyp2_1 + stat_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

5.4 Descriptive Effect Sizes Hyp 3/4

Spearman and Kendall correlations:

r_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)
GB CC CB
\(r(integrity, multiplism)\) -0.05 -0.13 0.05
\(\tau(integrity, multiplism)\) -0.04 -0.07 0

5.5 Hyp 4

data_scales_lables%>%
  mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
                               treat == "CB" ~ "Colored badges",
                               treat == "CC" ~ "Control condition",
                               T ~ "treat"))%>%
  ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
  geom_violin(adjust = 1.5, alpha = .5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  labs(title = "Graphical overview (Hyp 4)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  xlab("") +
  ylim(1,4) +
  hrbrthemes::theme_ipsum_ps()

fig4 <- data_scales_lables%>%
 mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
                              treat == "CB" ~ "Colored badges",
                              treat == "CC" ~ "Control condition",
                              T ~ "treat"))%>%
 ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
 geom_violin(adjust = 1.5, alpha = .5) +
 stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
 coord_flip() +
 xlab("") +
 ylim(1,4) +
 hrbrthemes::theme_ipsum_ps()

5.6 Descriptive Effect Sizes Hyp 4

A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)


A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "GB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)


A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CC")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)
GB CC
CC A = 0.49, d = -0.02
CB A = 0.46, d = -0.15 A = 0.46, d = -0.13

6 Inference statistics (Hyp 1)

6.1 Description of the variables

All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:

  • GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open Practices
  • CC means, that the participants of the study could not learn if or if not the authors of the study used Open Practices
  • CBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open Practices

Finally, the variable session identified the study participants.

If we look descriptively at these variables:

data_scales_wide%>%
  my_skim(.)
Data summary
Name Piped data
Number of rows 257
Number of columns 13
_______________________
Column type frequency:
character 1
numeric 12
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
session 0 1 64 64 0 257 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist skewn kurto maxabszvalue
Benevolence_CB 96 0.63 4.94 1.17 1.75 4.00 5.00 5.75 7 ▁▃▇▆▆ 0.03 -0.66 2.74
Benevolence_CC 87 0.66 4.96 1.22 1.25 4.00 5.00 6.00 7 ▁▂▇▇▇ -0.34 -0.04 3.05
Benevolence_GB 74 0.71 4.68 1.20 1.00 4.00 4.75 5.50 7 ▁▂▇▇▅ -0.06 -0.21 3.07
Experitse_CB 96 0.63 5.11 1.28 1.00 4.00 5.00 6.00 7 ▁▂▆▆▇ -0.33 -0.38 3.21
Experitse_CC 87 0.66 5.21 1.23 1.00 4.33 5.33 6.00 7 ▁▁▅▆▇ -0.63 0.40 3.43
Experitse_GB 74 0.71 4.99 1.31 1.00 4.00 5.00 6.00 7 ▁▂▇▇▇ -0.31 -0.37 3.04
Integrity_CB 96 0.63 4.96 1.15 1.25 4.00 5.00 5.75 7 ▁▂▆▇▅ -0.11 -0.21 3.24
Integrity_CC 87 0.66 4.97 1.19 1.00 4.00 5.00 6.00 7 ▁▁▇▇▆ -0.30 0.21 3.33
Integrity_GB 74 0.71 4.73 1.22 1.00 4.00 4.50 5.75 7 ▁▂▇▆▅ -0.11 -0.13 3.06
TSM_CB 96 0.63 2.60 0.58 1.00 2.25 2.50 3.00 4 ▁▃▇▃▂ 0.09 0.06 2.74
TSM_CC 87 0.66 2.67 0.57 1.00 2.25 2.75 3.00 4 ▁▂▇▃▂ 0.04 0.20 2.93
TSM_GB 74 0.71 2.68 0.55 1.00 2.25 2.75 3.00 4 ▁▂▇▅▂ 0.14 0.13 3.08

6.2 Investigating the missingness

6.2.1 Missingness per Variable

library(naniar)
## 
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
## 
##     n_complete
visdat::vis_miss(data_scales_wide%>%select(-session)) + 
  ggtitle("Missingness per Variable") + 
  theme(plot.margin = margin(1, 2, 1, 1, "cm"))

6.2.2 Marginal distributions Integrity_GB

library(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.3 Marginal distributions Integrity_CC

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.4 Marginal distributions Integrity_CB

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.5 Marginal distributions TSM_GB

ggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.6 Marginal distributions TSM_CC

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.7 Marginal distributions TSM_CB

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.8 Cohen’s d of missing/not-missing per variable

d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
  round(.,4)

knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")
Boxplot of Cohen’s d of missing/not-missing per variable
Integrity_CB Integrity_CC Integrity_GB TSM_CB TSM_CC TSM_GB
Integrity_CB NaN -0.1487 0.1183 NaN -0.0571 0.0617
Integrity_CC -0.0627 NaN -0.1183 -0.1243 NaN -0.0617
Integrity_GB 0.0627 0.1487 NaN 0.1243 0.0571 NaN
TSM_CB NaN -0.1487 0.1183 NaN -0.0571 0.0617
TSM_CC -0.0627 NaN -0.1183 -0.1243 NaN -0.0617
TSM_GB 0.0627 0.1487 NaN 0.1243 0.0571 NaN
boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")

6.3 Imputation

M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
            m = M,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

6.4 Check of first 10 imputation

out_first10 <-   mice(data = data_scales_wide%>%select(-session),
            m = 10,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

densityplot(out_first10)

plot(out_first10)

6.5 Parameter and BF estimation

# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3) 
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3) 


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance


# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))* 
  sum(diag(covbetween_hyp1 %*% 
             MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)


# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp1 <- bain(estimates_hyp1,
               "Integrity_GB<Integrity_CC<Integrity_CB;
                Integrity_GB=Integrity_CC=Integrity_CB;
                Integrity_GB<Integrity_CC=Integrity_CB",
                n = neff_hyp1, Sigma=covariance_hyp1,    
                group_parameters=3, joint_parameters = 0)
print(results_hyp1)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.516 0.161 3.199  5.543  0.135 0.130
## H2 0.476 0.237 2.012  2.012  0.085 0.082
## H3 4.011 0.217 18.473 18.473 0.780 0.748
## Hu                                 0.041
## 
## Hypotheses:
##   H1: Integrity_GB<Integrity_CC<Integrity_CB
##   H2: Integrity_GB=Integrity_CC=Integrity_CB
##   H3: Integrity_GB<Integrity_CC=Integrity_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
summary(results_hyp1)
##      Parameter        n Estimate       lb       ub
## 1 Integrity_GB 172.6037 4.746086 4.581740 4.910432
## 2 Integrity_CC 172.6037 4.958640 4.794640 5.122640
## 3 Integrity_CB 172.6037 4.965926 4.802789 5.129063
results_hyp1$BFmatrix
##           H1       H2        H3
## H1 1.0000000 1.589591 0.1731532
## H2 0.6290927 1.000000 0.1089294
## H3 5.7752334 9.180258 1.0000000

7 Inference statistics (Hyp 2)

7.1 Parameter estimation

path_mod <- "Integrity_GB ~ TSM_GB
             Integrity_CC ~ TSM_CC             
             Integrity_CB ~ TSM_CB"

# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 path_fit <- sem(path_mod, 
                 data = mice::complete(out, i), 
                 std.lv = T) # estimate the path coefficients 
 best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
                     filter(op == "~")%>%
                     pull(std.all) 
 covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
                  1/M * lavInspect(path_fit, 
                                   "vcov.std.all")[c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB"),
                                                   c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB")]
}

# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)
## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB 
##         -0.03         -0.13         -0.09

7.2 Visual path model

\usetikzlibrary{arrows,positioning}
\usetikzlibrary{calc}
\begin{tikzpicture}[auto,>=latex,align=center,
latent/.style={circle,draw, thick,inner sep=0pt,minimum size=10mm},
manifest/.style={rectangle,draw, thick,inner sep=2pt,minimum size=15mm},
manifestcov/.style={rectangle,draw, thick,inner sep=1pt,minimum size=8mm},
mean/.style={regular polygon,regular polygon sides=3,draw,thick,inner sep=0pt,minimum size=10mm},
paths/.style={->, thick, >=latex, font = \footnotesize\sffamily, fill = white},
variance/.style={<->, thin, >=latex, bend left=270, looseness=2.5, font = \footnotesize},
covariance/.style={<->, thin, >=latex, bend left=290, looseness=.5, font = \footnotesize},
dotsdots/.style={circle, draw, fill = black, minimum size = 1mm, maximum size = 1mm}
]

%% AVs
\node [manifest] (AV1) at (0, -2)       {$int_{gb}$};
\node [manifest] (AV2) [right = of AV1] {$int_{cc}$};
\node [manifest] (AV3) [right = of AV2] {$int_{cb}$};


%% UVs
\node [manifest] (UV1) at (0, 2)        {$tsm_{gb}$};    
\node [manifest] (UV2) [right =of UV1]  {$tsm_{cc}$};    
\node [manifest] (UV3) [right =of UV2]  {$tsm_{cb}$};    

%% Paths
\draw [paths] (UV1) to node [midway]{-.03} (AV1);
\draw [paths] (UV2) to node [midway]{-.13} (AV2);
\draw [paths] (UV3) to node [midway]{-.09} (AV3);

%\draw [covariance] (UV2) to node {} (UV1);
\end
{tikzpicture}

covbetween_hyp2 <- cov(best_hyp2) # between covariance matrix
covariance_hyp2 <- covwithin_hyp2 + (1+1/M)*covbetween_hyp2 # total variance

# Determine the effective and real sample sizes ######
samp_hyp2 <- nrow(data_scales_wide) # real sample size
nucom_hyp2 <- samp_hyp2-length(estimates_hyp2)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp2 <- (1+1/M)*(1/length(estimates_hyp2))* 
  sum(diag(covbetween_hyp2 %*% 
             MASS::ginv(covariance_hyp2))) # ... (43)
nuold_hyp2 <- (M-1)/(lam_hyp2^2) # ... (44)
nuobs_hyp2 <- (nucom_hyp2+1)/(nucom_hyp2+3)*nucom_hyp2*(1-lam_hyp2) # ... (46)
nu_hyp2 <- nuold_hyp2*nuobs_hyp2/(nuold_hyp2+nuobs_hyp2) # ... (47)
fracmis_hyp2 <- (nu_hyp2+1)/(nu_hyp2+3)*lam_hyp2 + 2/(nu_hyp2+3) # ... (48)
neff_hyp2 <- samp_hyp2-samp_hyp2*fracmis_hyp2 # = 114 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp2 <- list(covariance_hyp2)

7.3 Bayes Factor estimation (Hyp 2)

set.seed(6346)
results_hyp2 <- bain(estimates_hyp2, 
                     "Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
                     Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = F)
print(results_hyp2)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.570  0.163 3.495  6.804  0.094 0.092
## H2 22.961 0.685 33.497 33.497 0.906 0.882
## Hu                                  0.026
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
##   H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp2$BFmatrix
##          H1        H2
## H1 1.000000 0.1043274
## H2 9.585207 1.0000000

8 Inference statistics (Hyp 3)

set.seed(6346)
results_hyp3 <- bain(estimates_hyp2, 
                    "Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
                     (Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
                     Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = T)

print(results_hyp3)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.050  0.169 0.295  0.258  0.009 0.009
## H2 0.259  0.346 0.750  0.663  0.023 0.023
## H3 10.821 0.347 31.201 31.201 0.968 0.938
## Hu                                  0.030
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
##   H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
##   H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp3$BFmatrix
##           H1         H2          H3
## H1   1.00000  0.3938263 0.009469693
## H2   2.53919  1.0000000 0.024045352
## H3 105.60005 41.5880783 1.000000000

9 Inference statistics (Hyp 4)

# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance


# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))* 
  sum(diag(covbetween_hyp4 %*% 
             MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)


# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp4 <- bain(estimates_hyp4,
               "TSM_GB>TSM_CC>TSM_CB;
                TSM_GB=TSM_CC=TSM_CB;
                (TSM_GB,TSM_CC)>TSM_CB",
                n = neff_hyp4, Sigma=covariance_hyp4,    
                group_parameters=3, joint_parameters = 0)

print(results_hyp4)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.619  0.176 3.527  7.642  0.125 0.120
## H2 17.926 0.799 22.432 22.432 0.792 0.765
## H3 0.826  0.350 2.359  8.796  0.083 0.080
## Hu                                  0.034
## 
## Hypotheses:
##   H1: TSM_GB>TSM_CC>TSM_CB
##   H2: TSM_GB=TSM_CC=TSM_CB
##   H3: (TSM_GB,TSM_CC)>TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp4$BFmatrix
##           H1        H2       H3
## H1 1.0000000 0.1572512 1.495238
## H2 6.3592510 1.0000000 9.508594
## H3 0.6687899 0.1051680 1.000000